Non-empirical 'derivation' of B88 exchange functional 
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The B88 exchange energy density functional (created by Becke in 1988) is a crucial part of the most 
popular density functional in use today, B3LYP. B88 contains one empirical parameter which was 
fitted to Hartree-Fock exchange energies for the noble gas atoms. We show how local approximations 
to exchange become relatively exact under a very specific approach to the limit of large numbers, but 
that the usual gradient expansion does not. The leading corrections can be captured by generalized 
gradient approximations, producing a non-empirical derivation of the parameter in B88. 



I. INTRODUCTION 

Density functional theory (DFT) has become the 
method of choice for many electronic structure calcula- 
tions in quantum chemistry It balances the demands 
of accuracy and computation time, making it advanta- 
geous to other available methods. Although DFT is a for- 
mally rigorous theory, in practice it requires an approx- 
imation to the exchange-correlation (XC) energy, E'xc, 
as a functional of the density. In the late 1980's and 
early 1990's, approximations, such as B88, were devel- 
oped which proved accurate enough for their application 
in many areas of both chemistry and physics. The explo- 
sion in the use of DFT, driven by newly available com- 
putational power, could not have succeeded without the 
development of such functionals. 

Modern DFT began with the Hohenberg-Kohn theo- 
rem of 19640 and the Kohn-Sham scheme of 1965[1|. 
The first approximation to the XC energy was the local 
density approximation (LDA), which uses only the den- 
sity at a given point to determine the energy density at 
that point. This is exact for a uniform gas, since it is the 
known energy density of the uniform gas that is used to 
define LDA. The natural successor to LDA is a semi-local 
(or gradient-corrected) approximation which adds infor- 
mation about the derivative of the density at that point. 
In fact, in the same paper in which LDA is introduced, 
so too is the gradient expansion approximation (GEA) 
for XC. The coefficients of the GEA are determined by 
the energy of a slowly- varying gas 0, H, § . However it 
was found that the GEA often worsened LDA results 
and two decades passed before substantial improvements 
were made. 

Generalized gradient approximations (GGAs) effec- 
tively resum the gradient expansion, but using only | Vrij. 
The BBS functional^ is the most used GGA for exchange 
overall (as part of B3LYP[1, Q), but the most popular 
GGA in solid state applications is PBE[i3|. Neither re- 
duces to the GEA in the limit of small gradients. In 
this paper we explain the reason why this must be the 
case. Asymptotic expressions for the energy components 



as functionals of A^, the number of electrons, display 'un- 
reasonable accuracy' [ll| even for low N . In order to give 
good energies for finite systems, any approximate XC 
functional must have accurate coefficients in its large- iV 
expansion. LDA gives the dominant contribution, but 
the GEA does not yield an accurate leading correction 
for atoms. Popular GGAs such as B88 and PBE do get 
this correction right. 

In Ref. the underlying ideas behind this work were 
developed, however the reasoning was based upon scaling 
the density and not on the potential scaling discussed 
below. We refine these ideas and explicitly show how they 
can be used for functional development, and in particular 
we show how the parameter in B88 may be derived in a 
non-empirical manner. 



II. BACKGROUND THEORY 

In this section we review the theory behind asymp- 
totic expansions of the energy, including the expression 
for the exchange energy. A demonstration of the use- 
fulness of this asymptotic expansion for exchange is also 
given, where the constant in LDA is found without re- 
ferring to the electron gas. We also show the form of 
the gradient correction to LDA exchange, as well as dis- 
cussing the form of the GGAs, B88 and PBE. 



A. Asymptotic expansion in A*' 

Begin with any system (atom, molecule, cluster, or 
solid) containing N electrons. We then imagine changing 
the number of electrons to A^'. Since we usually begin 
from a neutral system, usually we consider only N' > N . 
Thus we define a scaling parameter = N' /N > 1. As 
we change the particle number, we simultaneously change 
the one-body potential Wcxt (r) in such a way as to retain 
overall charge neutrality, which means 



t(r) = C'/3^;,,t(C'/'r), N^CN. 
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We refer to this as charge-neutral (CN) scahng. For an 
isolated atom, Z — » C,Z under this scaling, so it remains 
neutral as the electron number grows. For molecules with 
nuclear positions and charges Zq,, — > C-^a and 
Rq C^^'''^R'a- In the special case of neutral atoms, 
the resulting series for the energy is well-known: 
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(2) 



where oq = 0.768745, ai = -1/2, and 03 = 0.269900 
pH . [Tsj . We say an approximation is large- asymptoti- 
cally exact to the p-th degree (AEp) if it recovers exactly 
the first p -I- 1 coefficients for a given quantity under the 
potential scaling of Eq. ([T]). Lieb and Simon[3, [Tsj 
showed that Thomas-Fermi (TF) theory becomes exact 
in the limit C ~^ 00 for all systems. TF is exact in a 
statistical sense, in that TF gives the correct first term 
of Eq. ([2]), but not the other terms. We say TF is AEO 
for the total energy. 

A similar expression exists for the exchange component 
of the energy alone: 
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where cq = 0.2208 = 9a2/ll and ci will be the main topic 
of this paper. In a similar fashion, Schwingcr demon- 
strated that the local approximation for exchange is AEO, 
and this coefficient is given exactl y by local exchange 
evaluated on the TF density [ll|, [H, [ij- However in or- 
der to give atomic exchange energies needed for chemical 
accuracy, any exchange approximation should be at least 
AEl. 

Now suppose we want to make a local approximation 
for Ex but know nothing about the uniform gas. Di- 
mensional analysis (coordinate scaling[l3l) tells us that 
it must be of the form: 



E]f'^[n]=AxI, 



1= d 



,4/3 



(r). 



(4) 



Requiring that this gives the leading contribution to Eq. 
([3]) then fixes the value of the constant A^. Using any 
(all-electron) Hartrec-Fock atomic code, such as were al- 
ready available in the 1960's[ll|, one calculates / for 
densities running down a particular column of the pe- 
riodic table and then deduces its dependence on Z^/^ . A 
modern alternative is to use the fully numerical OPMKS 
code[l^ using the OEP exact exchange functional to find 
densities for neutral atoms from Z = 1 to Z = 88. By fit- 
ting, one finds / = 0.2965^^/3 and hence Ax = -0.7446. 
This is remarkably close to the derived result of 



A. 



ran 1/3 



-0.7386. 



(5) 



Thus, without any recourse to the uniform electron gas, 
we have derived the correct local approximation to E^ [n] . 
This demonstrates that, via asymptotic exactness, the 
local approximation to exchange is a universal feature of 
all systems as ^ 00, when scaled appropriately. (In 
fact, Schwinger only proved this for atoms [llj. we know 
of no proof for arbitrary systems). 



B. Gradient expansions 

The next step up the ladder of increasingly sophis- 
ticated density- functional approximations [lOl is a semi- 
local approximation for E^ln], i.e., one which includes 
information about the gradient of the density. We intro- 
duce the dimensionless measure of the gradient: 
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|Vn(r) 
2fcp(r)7i(r) 
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where kp{r) = (37r^n(r))-'^/'^ is the local Fermi wavevec- 



tor. This is often written in terms of 



4/3 



which is simply proportional to s. Assuming smoothness 
in s and no preferred spatial direction, we know any sen- 
sible approximation depends only on s^. The gradient 
expansion is defined as the expansion of the energy as a 
functional of the density around the uniform limit. The 
leading correction for exchange is: 



E^^\n]^^,Jd'rs\v)e}f^ 



(n(r)). 



(7) 



where ex = A-^n'^^^ and ^ is a constant. Alternatively, 
we may write: 



4'^W = -/3 J < 



with 
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In a very slowly-varying electron gas, the gradient is very 
small, and the exchange energy will be accurately given 
by E^^^ + E^\ For such systems, the constant /i = 
lO/SlQ, so that (3 « 0.0024. 

The gradient expansion approximation (GEA) means 
applying this form to a finite system, using the value of 
ji from the slowly-varying gas. The GEA for exchange 
typically reduces the LDA error by about 50%. However 
it's counterpart for correlation worsens the LDA error, as 
its energy density is not even always negative. In many 
cases, GEA strongly overcorrects LDA leading to positive 
correlation energies and giving poor total energies 

A generalized gradient approximation (GGA) seeks to 
include the information contained in s(r) while improving 
on the success of LDA. The B88 exchange functional was 
designed to reduce to the GEA form when s is small, but 
also recover the correct —n{r)/2r decay of the exchange 
energy density for large r in atoms. Thus it interpolates 
between two known limits, and has the form: 



l-h6a;/3^«*^ sinh"^ [2i/3a;] ' 
(10) 

where AE^ denotes the correction to LDA. Thus the 
B88 functional^ contains one unknown parameter, /3^****. 




FIG. 1: The OEP exact exchange energies Exx for neutral 
atoms from Z = 1 to 88, divided by Z^^^ in order to pick 
out the leading term in its asymptotic series. The leading 
corrections are proportional to powers of ^"^/^ The values 
for the noble gas atoms are given as the circle symbols. 




FIG. 2: We add to Fig. [T]the resuhs for the LDA functional 
evaluated on the OEP densities (dashed-line) keeping the XX 
values (solid line). As in Fig. \T\ the Noble gas atom are 
highlighted with circle and square symbols for XX and LDA 
respectively. 



In 1988 Bcckc found this parameter by fitting to the 
Hartrcc-Fock exchange energies of the noble gases, find- 
ing a value of 0.0053. In fact, Beeke notes that this value 
is consistent with the observation of a high-Z asymptote 
for p. In Ref. [U, Becke calculates what value of (3 in 
Eq. ([5]) is required in order to give the HF exchange en- 
ergy for each atom in the first two rows of the periodic 
table along with the noble gas atoms. Thus, P is treated 
as a function of Z, and he observes that it converges for 
high-Z. Thanks to the previous section on asymptotic se- 
ries, we can now understand why this convergence occurs. 
Although the B88 form reduces to that of the gradient 
expansion for small gradients, the value for (3 is about 
twice as large as that predicted from the slowly-varying 
gas. 

Another common GGA for exchange is the Perdew- 
Burke-Ernzerhof (PBE) approximation [Toj . usually writ- 
ten in terms of an enhancement factor, Fx{s), to the LDA 
exchange energy density: 

<^^[n]=|dVFr"(.).^°^N (11) 

where 



and n = 0.2195 and k = 0.8040. This form for the en- 
hancement factor is chosen so that it reduces to LDA 
for s = and again recovers the form of the gradient 
expansion for small s. For large s it becomes a con- 
stant determined by the parameter k. Both k and fi 
are determined via satisfaction of various exact condi- 
tions. The value of fi was chosen to preserve the good 
linear response of LDA for the uniform electron gas un- 
der a weak perturbation [13, HI], while k is set by the 



Lieb-Oxford bound [IJl on the exchange-correlation en- 
ergy. (That condition is obviously violated by B88, while 
PBE does not accurately recover the X energy density in 
the tails of Coulombic systems). 



III. THEORY 

LDA yields the dominant term in either the asymp- 
totic charge-neutral expansion (C —>■ oo) or the gradient 
expansion for the slowly-varying electron gas, s — > 0. We 
next show that, contrary to popular myth, the impor- 
tant expansion is the charge-neutral expansion, not the 
gradient expansion. 

The charge-neutral expansion can be applied to any 
type of matter, be it molecule or extended solid. For 
any finite system, the density decays exponentially far 
from the nuclei. This is a key distinction between finite 
systems and bulk matter, treated with periodic boundary 
conditions. Bulk matter has no such regions. 

But, for slowly- varying gases, or more generally when 
there are no classical turning points at the Fermi surface, 
the charge-neutral scaling and the gradient expansion be- 
come identical, i.e., the gradient expansion for the slowly- 
varying gas is simply a special case of charge-neutral scal- 
ing. To sec that this is so, consider just the kinetic energy 
as a density-functional. Here the gradient expansion is 
known out to 6th order [25|. and eventually the integrated 
quantity itself diverges for atomic densities, due to the 
evanescent tail. But no such divergence occurs for ex- 
tended systems with finite density everywhere [25l l26t . 

Thus CN scaling applies to all systems, but only be- 
comes identical to the gradient expansion for slowly- 
varying bulk systems. For the dominant contribution, ef- 
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FIG. 3: We now find the next coefiicient in the exchange 
asymptotic series. To minimize the error due to shell structure 
oscillations, the LDA exchange energy is subtracted from the 
exact exchange for each atom. As both give the leading cor- 
rection, their difference will then have AcZ, Ac = Ci — c\^^, 
as the leading term in its asymptotic expansion, . The dashed 
line is the result of fitting to the noble gas atoms (circle sym- 
bols). 

fectively the local Fermi wavelength becomes short on the 
length scale on which the density is changing, so that the 
local approximation applies, and yields the exact answer 
for this term. Hence LDA reproduces the AEO terms, 
but GEA does not produce the leading corrections. All 
this has been amply demonstrated for simple Id model 
svstems[l6ll. and for the Kohn-Sham kinetic energy for 
atoms [231 • 

Here we apply the same reasoning for exchange. The 
local approximation becomes relatively exact as C cxd, 
but the gradient expansion docs not reproduce the lead- 
ing correction in the CN expansion. Below, we use the 
simple reasoning of Ref. [l^ to recover this leading cor- 
rection. We perform a much more extensive calculation 
of the asymptotic behavior, using methods developed in 
Ref. [13. We find, in agreement with Ref. [H, that the 
leading correction for atoms is about double that given 
by the gradient expansion, matching quite closely that of 
B88 and of PBE. Reversing this logic for B88, we show 
that B88 may be more or less derived non-empirically 
via the constraint that the approximation be AEl. If we 
enforce AEl exactly, we find a slightly different value for 
/?, and discuss the properties of the resulting functional, 
excogitated B88. 

IV. EXTRACTING ASYMPTOTIC 
COEFFICIENTS 

Under the potential scaling of Eq. ([1]), any approxima- 
tion for the exchange energy that reduces to LDA in the 
uniform limit has an expansion in N like Eq. ([3]), with 




20 40 60 

z 



FIG. 4: In order to see that LDA does not significantly con- 
tribute to the higher orders of the exchange asymptotic series, 
we plot the difference between LDA and the leading term 
cqZ^^"^ (dashed line), as a function of Z. The exact exchange 
value is also shown (solid line). 



the same value for cq. However the coefficient ci depends 
on the particular approximation. Below we explain the 
procedure used to extract these coefficients. 

As mentioned in the previous section, the OPMKSp^ 
electronic structure code is a fully numerical electronic 
structure code that has the ability to perform optimized 
effective potential (OEP) calculations. We evaluate the 
various approximations using atomic densities found with 
the OEP exact exchange (XX) method. The densities 
found using this method will be extremely close to the 
exact densities despite the fact that correlation is miss- 
ing. Moreover the effect of correlation will contribute at 
higher orders in the asymptotic expansions of the energy 
than those we are interested in. Thus EXX calculations 
arc in principle sufficient for extracting the coefficient we 
seek. 

In Fig. [2 we plot E^^/Z^/^ vs Z-^/^ where E^^ is 
the exchange energy from the exact exchange calculation. 
Since the leading term in the asymptotic expansion of Eq. 
^ is Z^^^, this procedure picks out the cq coefficient as 
a constant while all other terms are functions of Z^^^"^. 
One can see that the curve in Fig. [T] is heading towards 
the exact value of Cq = —0.2208, but it is difhcult to 
extract higher coefficients due to oscillation of the curve 
due to the shell structure. 

To overcome this difficulty, in Fig[2] we add the LDA 
curve to Fig[TJ It can be seen that it too recovers the ex- 
act Co coefficient, but also clearly differs in higher orders 
in the asymptotic expansion. More usefully, we see that 
the LDA curve mimics the oscillations shown by exact 
exchange, so subtracting LDA from EXX will minimize 
this effect and make the extraction of asymptotic coeffi- 
cients more accurate. 

In Fig. [3 we plot (Sxx - E^^^)/Z vs Z^^/^ and find 
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FIG. 5: We use the same procedure as in Fig. [3] to find the 
Ac coefficient for the gradient correction to LDA, E^^^ [n] , as 
defined in Eq. (O. The dashed line was fitted to the noble 
gas atoms (circle symbols). 



FIG. 6: The percentage error of the approximate asymptotic 
series given in Eq. (|13|l is plotted as a function of Z. The 
error is remarkably low and demonstrates the power of these 
asymptotic series. 



that it behaves elose to Unearly. There appears to be no 

^4/3 

term in Ey^. Such a contribution was argued not 
to exist in Ref. [13, but this was based on CN scahng of 
the density and studying the behavior of the terms in the 
gradient expansion. That reasoning is insufficient, as the 
expansion should be performed in terms of the potential, 
as described in Sec. Ill Al But since the Scott correction 
(the contribution to the total energy) comes from the 
core region, there is no reason to expect an analogous 
contribution for exchange. In order to show this and also 
to precisely determine the ci coefficient, one should use 
the techniques developed by Schwinger for deriving the 
Scott correction to the total energy llj , but apply them 
to exchange. 

To further reduce the remaining uncertainty due to 
shell structure oscillations, we choose simply to use the 
noble gas atoms (excluding helium) for our fit. The 
strongest deviations from linearity come from the tran- 
sition metals and lanthanides and actinides. We fit the 
difference {E^y^- E]l^^) / Z with a straight line in Z''^/^, 
and extrapolate to Z — > oo, finding Ac = —0.2240, where 
Ac = ci - c\^^, and 

-Bxx ~ E^^^ - 0.2240 Z + 0.2467 Z^/^. (13) 

The coefficient of the last term is given by the slope of 
the dashed line in Fig. [31 although the meaning of this 
term is unclear in the presence of such strong oscillating 
contributions. 

If instead we used the alkaline earth atoms (exclud- 
ing beryllium), we find an almost identical value, Ac = 
—0.2236. If we use all elements with Z > 10, we find a 
similar value, Ac = —0.2164. If all elements from Z = 1 
to 88 are used, we find Ac = -0.1982. In Ref. [H, Ac 
was found using noble gas atoms, except with the he- 
lium value included, and that method gave a value of 



Ac = —0.1978. In our analysis, atoms with Z < 10 are 
not used as they are not necessarily dominated by the 
asymptotic series. 

Since LDA displays the shell oscillations that pre- 
vented us from fitting EXX directly, the value of the LDA 
ci coefficient cannot be found exactly. But we estimate 
> c\^^ > —0.04, i.e., at least five times smaller in 
magnitude than Ac. In Fig[4l we show E^ — cq Z^^^ as a 
function of Z for both the exact values and within LDA, 
demonstrating that the linear contribution comes almost 
entirely from the beyond-LDA terms. 

Finally we determine Ac for GEA. In Fig. [51 we 
plot (E^EA _ £;LDA)/^ 2-1/3 in order to find Ac = 
^GEA _ ^LDA^ finding Ac = -0.1062. This plot is much 
closer to linear than the previous one. The leading cor- 
rections to LDA in the asymptotic expansion produce 
corrections to the shell structure beyond those captured 
by LDA evaluated on the exact density [IB, [131 • Although 
the smooth contribution can be partially captured by 
GEA, there is almost no correction to the shell struc- 
ture. Just as for the kinetic energv[27|. GEA yields a 
correction to the smooth part that is about half of the 
accurate value. 



To understand the power of these asymptotic expan- 
sions, we add the corrections of Eq. to the LDA 
energies, and in Fig. [6] plot the percentage error relative 
to exact exchange, as a function of Z. For all but the 
second row of the periodic table, the resulting error is 
below 0.5% in magnitude, and typically of order 0.2%. 
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FIG. 7: The percentage error for LDA, GEA and the mod- 
ified GEA (MGEA) exchange functional for Z > 10. The 
coefficient of i^x is muhipUed by 2.109 in order to make the 
MGEA AEl asymptoticaUy exact. 



TABLE I: Ac = ci - ci values for several diff'erent func- 
tionals. 





Exx 


LDA 


GEA 


B88 


PBE 


Ac 


-0.2240 




-0.1062 


-0.2216 


-0.1946 



V. GENERALIZED GRADIENT 
APPROXIMATIONS 

Generalized gradient expansions were designed to im- 
prove energetics over LDA for electron systems of inter- 
est and relevance. Early versions, such as PW9l[28l[29t. 
were tortured into reducing to the gradient expansion 
when the density is slowly varying. But this was later 
given up, in both B88 and PBE exchange, which both re- 
duce to the gradient expansion form for slow variations, 
but with coefficients much larger than that of the gradi- 
ent expansion. 

Our analysis explains why this must be so. Regardless 
of its derivation, any modern GGA for exchange is tested 
against the neutral atoms. Any approximation that can- 
not recover the right c\ will be generally inaccurate for 
these energies, and so discarded. Thus any that become 
popular have already passed this test. 

In Table HI we give the results for Ac for several dif- 
ferent functionals. The same methodology was used in 
all extractions. Both popular GGA's recover (at least 
approximately) the accurate value. B88, designed specif- 
ically for molecular systems, is very close to the accurate 
value. PBE exchange is less so, but is also designed to 
bridge molecular and solid-state systems. The PBE value 
is between that of GEA and B88, but much closer to the 
latter than the former. Taking advanta ge o f this insight, 
a new variation on PBE, called PBEsol[30|, restores the 
original gradient expansion, thereby worsening atomiza- 



FIG. 8: We add to Fig. [7] the resufis for B88, PBE, and the 
excogitated B88 functional, all evaluated on the OEP exact 
exchange densities. 



tion energies (total energy differences), but improving 
many lattice constants of solids over PBE and LDA. 



A. Deriving the /3 in B88 

The exchange energies found using GGAs such as B88 
or PBE are generally dominated by their gradient expan- 
sion components for most chemically relevant densities. 
Thus, to make B88 AEl, it is sufficient to impose this ex- 
act condition on just the E'x functional form. Since 
both B88 and the GEA are built on top of LDA, we can 
simple look at the Ac values calculated above. If we set 
/i = 2.109/LiAK in Eq. (I7|), we multiply the c\ coefficient 
of GEA by a factor of 2.109, making it AEl. In Fig. [3 
we name this functional MGEA for modified GEA and 
plot its percentage errors. The values for LDA and GEA 
are also shown. It can be seen that modifying GEA to 
be AEl has greatly reduced the error. 

We now require that the B88 functional form, Eq. (|10|) . 
reduce to this MGEA for small values of x. Using Eq. 
([5]), this corresponds to using a value of /? = 0.0050. Thus 
we have derived an excogitated B88 that is free of any 
empirical parameters. The actual value used in B88 is 
^B88 „ 0.0053 (for spin-polarized systems this becomes 
0.0042, which is the value given in Ref. (0)), so the val- 
ues are very close. This is not surprising as fitting to 
Hartree-Fock exchange energies is an approximate way 
of demanding asymptotic exactness. Interestingly, the 
value quoted as the high-Z asymptote in Ref. and 
found using Ref. [2l| is essentially the same as our value, 
but was evidently rejected in favor of a better fit. In Fig. 
m we plot the percentage errors for B88, PBE, and the 
excogitated B88. As is typical for empirically-fitted func- 
tionals, B88 performs very well for systems close to the 
data set used in the fitted procedure. Although the er- 
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ror for PBE is higher than both B88 and the excogitated 
B88, it is systematic in its overestimation. As noted, 
PBE was designed to perform reasonably well for a wide 
range of systems, so again its behavior is not surprising. 
On this data set, the excogitated B88 was never going to 
do better than B88, although it remains to be seen how 
it performs for more complicated systems. 

VI. CONCLUSION 

We have carefully and systematically extracted the 
leading large-Z correction to the exchange energy of 
atoms. Our results differ slightly from those of Ref. 
IT^ but yield the same qualitative conclusion, i.e., that 
the gradient expansion yields an error of a factor of 2 or 



more for this coefficient. We have clarified some of the 
reasoning, and applied it more generally to any atom, 
molecule, or cluster. By looking in detail at the ex- 
change energy asymptotic series for neutral atoms, we 
have demonstrated the power of using such series for 
functional development. Requiring that the small gradi- 
ent expansion of B88 capture the two leading coefficients 
of the asymptotic expansion is a method by which the 
unknown coefficient (3 can be found. This gives a co- 
efficient very close to the one actually used in B88 and 
thus is an ex post facto 'derivation' of B88. Inserting our 
most accurate estimate for /3 into the B88 form yields an 
excogitated B88. 

Wc thank Ebcrhard Engel for the use of his atomic 
OPMKS code, and NSF CHE-0809859 . 
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